home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / program / nrpas13.zip / BRENT.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  79 lines
  1. FUNCTION brent(ax,bx,cx,tol: real; VAR xmin: real): real;
  2. (* Programs using routine BRENT must supply an external function
  3. func(x:real):real whose minimum is to be found. *)
  4. LABEL 1,2,3;
  5. CONST
  6.    itmax=100;
  7.    cgold=0.3819660;
  8.    zeps=1.0e-10;
  9. VAR
  10.    a,b,d,e,etemp: real;
  11.    fu,fv,fw,fx: real;
  12.    iter: integer;
  13.    p,q,r,tol1,tol2: real;
  14.    u,v,w,x,xm: real;
  15. FUNCTION sign(a,b: real): real;
  16.    BEGIN
  17.       IF (b > 0.0) THEN sign := abs(a) ELSE sign := -abs(a)
  18.    END;
  19. BEGIN
  20.    IF ax < cx THEN a := ax ELSE a := cx;
  21.    IF ax > cx THEN b := ax ELSE b := cx;
  22.    v := bx;
  23.    w := v;
  24.    x := v;
  25.    e := 0.0;
  26.    fx := func(x);
  27.    fv := fx;
  28.    fw := fx;
  29.    FOR iter := 1 TO itmax DO BEGIN
  30.       xm := 0.5*(a+b);
  31.       tol1 := tol*abs(x)+zeps;
  32.       tol2 := 2.0*tol1;
  33.       IF (abs(x-xm) <= (tol2-0.5*(b-a))) THEN GOTO 3;
  34.       IF (abs(e) > tol1) THEN BEGIN
  35.          r := (x-w)*(fx-fv);
  36.          q := (x-v)*(fx-fw);
  37.          p := (x-v)*q-(x-w)*r;
  38.          q := 2.0*(q-r);
  39.          IF (q > 0.0) THEN  p := -p;
  40.          q := abs(q);
  41.          etemp := e;
  42.          e := d;
  43.          IF((abs(p) >= abs(0.5*q*etemp)) OR (p <= q*(a-x))
  44.             OR (p >= q*(b-x))) THEN GOTO 1;
  45.          d := p/q;
  46.          u := x+d;
  47.          IF (((u-a) < tol2) OR ((b-u) < tol2)) THEN d := sign(tol1,xm-x);
  48.          GOTO 2
  49.       END;
  50. 1:      IF (x >= xm)  THEN e := a-x ELSE e := b-x;
  51.       d := cgold*e;
  52. 2:      IF (abs(d) >= tol1)  THEN u := x+d ELSE u := x+sign(tol1,d);
  53.       fu := func(u);
  54.       IF (fu <= fx)  THEN BEGIN
  55.          IF (u >= x)  THEN a := x ELSE b := x;
  56.          v := w;
  57.          fv := fw;
  58.          w := x;
  59.          fw := fx;
  60.          x := u;
  61.          fx := fu
  62.       END ELSE BEGIN
  63.          IF (u < x)  THEN a := u ELSE b := u;
  64.          IF ((fu <= fw) OR (w = x))  THEN BEGIN
  65.             v := w;
  66.             fv := fw;
  67.             w := u;
  68.             fw := fu
  69.          END ELSE IF ((fu <= fv) OR (v = x) OR (v = 2)) THEN BEGIN
  70.             v := u;
  71.             fv := fu
  72.          END
  73.       END
  74.    END;
  75.    writeln('pause in routine BRENT - too many iterations');
  76. 3:   xmin := x;
  77.    brent := fx
  78. END;
  79.